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USING  LIDAR  TOMOGRAPHY  TO  CHARACTERIZE 
IN  A  TURBULENT  MEDIA 

1 .  Introduction 

Atmospheric  lidar  returns  depend  on  both  backscatter( P )  and 
extinction ( o )‘: 

P(r)=.^^,-^Uxp[-2  f  a(r')dr'],  (1) 

r'-o 

where  P  is  the  power  return,  r  is  range,  and  A*  contains  system 
constants.  Exact  determination  of  extinction  and  backscatter  from 
a  single  return  requires  a  knowledge  of  the  range  dependence  of 
either  extinction  or  backscatter  or  the  relationship  between  these 
variables.  However,  for  multiple  returns^,  tomography  methods  can 
be  used  to  independently  determine  extinction  and  backscatter.  In 
one  approach,  the  logarithmic  form  of  the  lidar  equation  separates 
the  two  variables: 

Sir)  «ln tP(r)  r^/X]  »ln[p  (r)  ]  -2t  (r) ,  (2) 

where  t  is  the  optical  depth  [T(r)«  f  o(r')dr')  and  the  lidar  data 

is  both  range  corrected  and  calibrated.  From  an  aircraft,  the  lidar 
can  scan  45*  forward  of  vertical,  straight  down,  and  45*  aft  of 
vertical,  see  figure  l: 

^45* 

5o.  “Ba-Tj  (3) 

where  B  and  T  are  respectively  the  logarithm  backscatter  [ln[0]  ] 


Figure  1.  An  aircraft  flying  from  left  to  right  makes  three 
tomographic  lidar  observations  of  the  same  volume. 


Manutcf^  appiovcd  Febniuy  23,  1994. 
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Controlled  scanning  allows  the  backscatter  at  the  one  point  to  be 
observed  several  times.  In  a  homogeneous  atmosphere,  the  system  of 
equations  becomes  over-determined: 

T 

^  * _ *2 

^  ^  COS(45“)  (4) 

and  the  backscatter  and  extinction  can  be  determined  from  a  set  of 
equations  that  use  a  least  square  fitting  process  to  minimize  the 
errors . 

While  equations(3)  and  (4)  are  overly  simplistic,  they 
illuminate  the  advantage  —  over-determined  equations  —  and 
limitation  of  lidar  tomography  —  the  need  for  homogeneity.  Spatial 
homogeneity  is  needed  since  the  scattering  must  be  almost 
homogeneous  inside  the  individual  grid  elements.  If  tomography  can 
use  regions  whose  scale  sizes  are  an  order  of  magnitude  smaller 
than  the  size  of  the  dominant  atmospheric  structures,  the 
approximation  can  be  quite  reasonable.  Temporal  homogeneity  is  also 
needed,  ie.  the  atmosphere  must  not  dramatically  change  between  the 
measurements.  Simply  stated,  the  time  between  the  Ildar 
measurements  must  be  substantially  faster  than  the  'half  life'  of 
the  dominant  atmospheric  structures. 

This  report  contains  an  initial  discussion  of  the  tomographic 
method  and  will  not  provide  a  detailed  analysis  of  equipment  and/or 
atmospheric  conditions  required  for  this  method  to  work.  While 
possible  instrumental  configurations  include:  mobile  lidars, 
multiple  fixed  lidars,  and  even  single  lidars,  only  a  single 
scanning  lidar  (either  moving  or  stationary)  is  modelled.  The  model 
used  to  test  the  inversion  accuracies  are  selected  for  simplicity 
and  ability  to  test  a  wide  range  of  conditions;  however,  the  model 
is  not  based  on  any  specific  atmospheric  model. 

2.  Basic  Technique  and  Bquaticms 

In  the  specific  tonography  method  discussed  by  this  report,  a 
grid  divides  the  atmosphere  into  equal  regions .  Bach  region  is 
assumed  to  have  a  constant  backscatter  and  extinction.  A  discrete 
form  of  the  lidar  equation  is  used: 

n.m 

i'.P-l.l 

where  i  and  j  represent  spatial  coordinates,  k  is  a  unique  number 
assigned  to  each  return,  n  and  m  are  the  maximum  size  for  the 
coordinates,  and  is  a  geometrical  factor.  For  extinction,  j 
is  determined  by  the  distemce  through  the  grid  element  that  the  kth 
profile  follows.  For  backscatter,  A^^  is  either  one  (if  this 
return  is  from  this  region)  or  zero.  A  series  of  measurements  of 
the  same  grid  are  made  from  different  locations,  systems  of 
equations  are  formed,  and  the  equations  are  coiabined  into  a  matrix 


relationship: 


where  the  matrix  equation  is  5  >  A  fi ,  K  is  the  total  number  of 
returns,  and  N  is  the  total  number  of  grid  elements  [n«/n]. 

Table  1  shows  the  grid  that  is  used  for  the  tomography 
calculations.  Each  square  is  arbitrarily  scaled  to  one  on  a  side. 
The  center  of  each  square  is  set  to  be  an  integer  and  the  edge  then 
becomes  an  integer  plus  (or  minus)  a  half.  The  wide  line  on  the 
left  hand  side,  which  represents  the  lidar  location,  has  a  j 
coordinate  of  0.5  .  The  first  set  of  observations  (or  scan) 
originate  at  the  top  [i=0.5]  and  the  other  sets  of  observations  are 
made  at  1/2  step  intervals  [ j«*0.5,l.0,1.5,  etc.].  Backscatter 
measurement  are  made  at  center  of  each  square  [  ie .  ( i»l .  0 ,  .  0 ) , 
( i*l . 0 , .0),....  ( i*2 . 0 , j*l .0),...]. 


backscattBring  and  the  bottom  number  is  the  matrix  number  for  the 
extinction. 
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The  coefficients  for  the  first  two  returns  are: 

[for  a  lidar  located  at  (i=0. 5, j=0. 5)  and  a  return  from 

(i=1.0, j=1.0)] 

(7) 

[In  the  above  case,  the  coefficients  ^  and  Aj  26  respectively  l.o 

and  All  other  coefficients  for  this  return  are  zero  (ie.A^  ^-0 

where  1  or  26  .  ] 

[for  a  lidar  located  at  (i=0.5, j=0.5)  and  a  return  from 

(i=1.0, j=2.0)] 

(8) 

[In  the  above  case,  the  coefficients  and  Aj.ji  are 

respectively  1.0,  ,  and  ] . 

The  scan  angles  are  restricted  to  less  than  or  equal  to  45°  [ie. 
the  colxunn  change  ( j)  is  less  than  or  equal  to  the  row  change  (1)]. 
This  limitation  is  not  necessary  and  represents  a  convenient,  but 
arbitrary  restriction.  For  the  grid  shown  in  ted>le  1  and,  for  the 
scheme  discussed  above,  195  returns  are  possible.  Since  the 
equations  have  50  unJcnotms,  the  equation  system  Is  over-determined 
and  the  least  squares  determination  of  the  atmospheric  parameters 
(backscatter  and  optical  depth)  becomes*: 

B*-a’*A-‘A=')^  (9) 

where  fi*  denotes  the  least  square  estimate,  A'  Is  the  transpose  of 
X,  and  A~^  Is  the  Inverse  of  A. 

The  added  Information  In  the  over-determined  equations  can  be 
used  In  several  ways.  Using  the  estimated  atmospheric  parameters, 
estimated  return  values  can  be  derived  and  an  error  [ e  ]  determined 
from  the  difference  between  the  actual  [5j(]  and  estimated  retxirn 

values 

e-E(5V-9*)^  (10) 

Then  this  difference  can  be  used  to  normalize  the  original 
geometrical  terms: 

where  the  weight  function  must  be  normalized  [  i  W^i-l  ] ,  and  the 

'  i'«i 

tomography  equations  recalculated  using  reduced  Input  errors. 
Another  use  of  the  extra  degrees  of  freedom  could  be  to  account  for 


temporal  changes  in  the  scattering  media.  For  example,  the 
resulting  terms  might  account  for  changes  in  optical  properties 
when  a  constant  dry  aerosol  size  distribution  is  modified  by 
changing  humidity. 

Computationally,  the  inversion  of  the  matrix  with  geometrical 
terms  requires  the  n‘  multiplications  [where  n  denotes  the  size  of 
measurement  grid]  and  can  not  be  easily  done  on  small  computers 
except  for  small  grids*.  However,  if  the  appropriate  inverse  of 
the  geometrical  matrix  [A]  is  calculated  ahead  of  time  on  a  large 
computer  and  transferred  to  a  small  computer,  the  number  of 
multiplications  is  reduced  to  n*.  Then  scattering  coefficients  can 
be  derived  on  the  PC  computer  from  equation  (9).  While  only  two 
dimensional  grids  are  considered  in  this  report,  a  three 
dimensional  grid  could  be  analyzed  with  the  same  techniques  (and 
such  measurements  are  well  within  the  capabilities  of  a  state-of- 
the-art  lidar) . 

3.  Kand<»  Errors  and  Error  Reduction 

The  tomography  analysis  technique  outlined  in  section  2  is  a 
'least  squares'  process  and  the  errors  in  the  derived  parameters 
matrix  are  linearly  related  to  the  magnitude  of  errors  in  the 
measurement  matrix: 

(12) 

«rt:ere  §*  and  S*  are  the  errors  in  the  measurement  and  derived 
matrices  respectively. 

To  test  the  sensitivity  of  this  inversion  algorithm,  a 
thousand  backscatter  and  extinction  matrices  were  randomly  selected 
with  no  relationship  between  backscatter  and  extinction.  The 
logarithm  of  the  backscatter[ values  are  selected  from  a 
Gaussian  distribution  with  an  average  of  zero  and  a  rms  variation 
of  one.  The  extinction  [T^],  which  corresponds  to  the  optical  depth 
of  each  square,  are  also  selected  from  a  Gaussian  distribution  with 
an  average  of  0.1/grid  element  and  a  variance  that  is  10%  of  the 
mean  value.  Using  the  selected  backscatter  and  extinction,  the 
simulated  return  values  [5]  are  derived  from  the  original 
geometrical  matrix  [if]  with  a  random  error  [and  with  a  variance  of 
5*]  added.  Table  2  shows  the  errors  in  the  extinction  and 
backscatter  terms  at  various  grid  locations: 


where  Eg  is  the  resulting  error,  J  is  the  total  number  of 
observations,  <md  the  original  input  for  atmospheric  parameters  and 
the  output  values  are  and  respectively.  Except  for  the  last 
column  [  j«5  ] ,  the  backscatter  errors  are  moderately  larger  than  the 


input  errors  with  only  a  small  variation  between  the  various  rows 
and  columns.  The  extinction  errors  are  larger  than  the  backscatter 
errors  and  increase  in  range.  Since  few  returns  from  the  last 
column  are  used,  the  backscatter  and  extinction  errors  are 
substantially  larger  for  this  column. 


Table  2:  Output  errors  resulting  from  uniform  input  errors 


top  and  bottom  line  of  each  box.  Output  errors  are  normalized  by 
input  errors. 

Table  3  shows  analysis  of  random  data  when  the  'return' 
[i«0.5, j»3.0]  from  the  center  square  [i»3.0, j*3.0]  has  random 
errors  added.  While  the  derived  values  from  the  grid  with  the  error 
are  influenced,  the  least  squares  process  spreads  the  errors  to 
other  squares  as  well.  The  last  column  again  has  the  largest  errors 
and,  in  the  case  of  extinction,  the  errors  in  the  center  square  are 
not  even  sxibstantially  larger  than  the  other  errors  in  the  center 
column. 

Using  the  derived  atmospheric  parameters,  new  'returns'  can  be 
calculated  from  the  geometrical  matrix  [A]  and  compared  with  the 
original  'returns'.  For  example,  in  the  case  discussed  above,  the 
difference  between  derived  and  original  returns  was  calculated.  The 
recalculated  return  (from  the  one  which  errors  were  added)  has  the 
largest  error,  which  reduces  the  original  error  by  80%.  Other 
recalculated  returns  from  the  sense  square  [l-0,j=:2.5  to 
i-3.5,j-3.5;  i-0,j-3.5  to  i»3.5, j«3.5;  i-0,j-4.5  to  i-3.5,j=3.5; 
i«0,j«4.0  tp  i«3.5,j-3.5]  increase  original  error  by  10%. 


6 


Table  3:  Output  errors  resulting  froa  an  error 
to  the  center  square 


Range 

[j] 

1 

a 

2 

3 

4 

5 

B 

m 

ra 

2.2 

3.8 

1.2 

10.2 

5.1 

2.2 

B 

Distance/ 

2 

m 

2.4 

4.3 

0.4 

9.3 

6.5 

22.8 

a 

Time  [ i ] 

3 

0.2 

2.3 

3.4 

4.5 

1.2 

10.9 

10.9 

8.0 

3 

B 

m 

m 

2.4 

4.3 

0.4 

9.3 

6.5 

22.8 

B 

iJ 

ra 

2.2 

3.8 

1.2 

10.2 

5.1 

2.2 

5 

2 

3 

4 

5 

top  and  bottom  line  of  each  box.  Output  errors  are  normalized  by 
input  errors  and  scaled  by  10'^. 


4.  Errors  introduced  by  a  changing  scattering  aedia 

Scattering  media  can  be  altered  by  redistribution  of 
scattering  structvures  and  changes  in  scattering  media  itself.  In 
this  section,  changes  caused  by  turbulence  to  scattering  structures 
are  simulated,  while  the  average  scattering  conditions  are  left 
unchanged.  Constant  movement  of  the  media  is  not  simulated  since  it 
is  assumed  that  movement  can  be  removed  during  the  initial  analysis 
of  the  lidar  data. 

The  spectrum  [F^]  of  scattering  structures  is  assumed  be 
Gaussian: 


F^(«l,«j)=n - -i-2 - , 


Where  u  is  the  wave  number,  the  numerical  subscript  is  the  grid 
component,  and  is  the  rms  width  of  the  spectrum.  Each  phase  of 
each  Fourier  component  is  selected  randomly  and  has  a  random  phase 
velocity: 

F(a>i,a>3,  t)  *Fo(»i,«2)exp{2wi  [riJdCWj,®,)  +(Vi+V2)  t]}, 

where  F  is  the  complex  magnitude  the  spectral  components,  t  is 
time,  phase  velocities  have  a  gaussian  distribution  e3q> [- (v^/o,.) ^1  , 
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and  the  wave  numbers  are  randomly  [  rnd((i)j^,  w^)  ]  selected.  The 
spectrum  with  appropriate  phases  is  calculated  on  a  16  by  16  grid 
and  Fourier  transformed  with  a  two  dimensional  FFT: 

t)  =£^j&exp[-2iiio)iXi/n]  exp  [-2nifc)2Xj/n]  F((Oi,  Wj,  t) , 

where  f  is  magnitude  of  the  scattering  in  the  real  domain  and  n  is 
the  maximum  grid  size  [16].  The  final  scattering  media  is  a  5  by  5 
sub^element  selected  from  the  larger  grid.  This  process  is  done 
twice  to  select  matrices  for  backscatter  and  matrices  for 
extinction.  Figure  2  shows  the  scattering  media  in  one  direction 
where  the  change  of  the  aerosols  structures  in  time  can  be  seen. 

Using  tomography,  backscatter  and  extinction  errors  are 
calculated  for  a  hundred  cases  each  with  eleven  time  steps.  These 
errors  (relative  to  the  original  errors)  were  dependant  on  the 
ratio  of  backscatter  to  extinction,  the  dominant  structure  size, 
and  the  rms  wind  variation  of  the  relative  structures.  Figure  3 
shows  the  backscatter  and  extinction  errors  for  the  different 
ratios  of  these  parameters.  The  smallest  errors  occur  when  the 
ratio  is  approximately  one.  When  the  extinction  values  are  small, 
the  error  in  the  backscatter  values  are  small  even  when  extinction 
errors  are  large.  Since  the  extinction  is  dependant  on  both 
backscatter  and  extinction,  this  parameter  is  strongly  dependant  on 
both  and  the  error  is  a  minimum  when  the  two  errors  are  comparable 
in  size. 


c 
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Figure  2:  The  changes  in 
a  typical  simulated 
scattering  media  are 
shown.  Each  line  shows 
the  scattering  in  the 
one  direction  through 
the  grid  at  one  time. 
The  different  fourier 
components  have 
different  phase 
velocities  which  are 
selected  from  a  Gaussian 
distribution  where  the 
rms  variation  equals  one 
grid  square  per  each 
time  element. 


8 


Figure  3 :  The  error 

scaled  for  random 
movement  and  structure 

size  is  plotted 

on  the  ordinate  for  the 
different  ratios  of 
backscatter  (natural 
logarithm)  and 
extinction.  The 
extinction  and 
backscatter  errors  are 
plotted  as  solid  and 

dashed  lines 
respectively.  Note,  the 
smallest  error  occurs 
when  the  ratio  is 

approximately  one. 


The  relative  errors  are  influenced  by  the  random  movement  of 
and  the  size  of  the  structures.  When  the  structures  move  rapidly, 
the  scattering  media  changes  quickly  and  the  variation  of  the 
scattering  media  over  the  measurement  cycle  becomes  larger.  Since 
the  smaller  structiires  move  a  larger  fraction  of  their  scale  size 
in  shorter  time  intervals,  the  relative  errors  are  also  dependant 
on  structure  sizes.  When  the  backscatter  to  extinction  ratio  is 
one,  the  relative  errors  are: 

P.rrl,^J2_ 

o„rl 

where  and  are  the  average  backscatter  and  extinction 

errors  divided  by  the  average  parameter  value,  o„  is  the  rms 
variation  of  speed  [in  units  of  one  grid  dimension  divide  by  one 
time  step],  and  is  the  rms  scale  size  of  aerosol  structures  [in 
units  of  one  grid  dimension]. 

If  each  grid  element  is  assumed  to  be  10  m  on  a  side  and 
random  movement  of  structxires  is  assumed  to  be  1  m/s,  the  error  for 
various  Ildar  configurations  can  be  calculated.  For  an  aircraft 
moving  at  100  m/s  over  a  media  which  has  a  dominant  structure  size 
*  of  1  km,  the  relative  error  becomes  15%,  which  is  acceptable  for 

many  atmospheric  measurements.  For  a  stationary  Ildar  when  the  wind 

speed  is  lOm/s,  the  error  is  150%  which  is  not  acceptable. 

«  ( 
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5.  Sunary  and  Conclusions 

Tomography  can  be  used  to  simultaneously  measure  both 
backscatter  and  extinction  measurements.  However,  the  method  is 
limited.  The  scattering  properties  need  to  be  almost  homogenous 
inside  the  tomography  grid  elements  and  must  not  change 
significantly  during  the  measurement  cycle.  Computationally  the 
most  multiplications  required  is  n* .  With  a  pre^calculated  inverted 
geometrical  matrix  [A],  the  solution  of  individual  terms  requires 
n*  multiplications. 

Since  the  equations  are  over-determined,  scattering  properties 
can  be  derived  from  a  least  squares  fitting  process.  Errors  are 
linearly  related  to  the  matrix  relationship,  v?hich  includes  linear 
extinction  and  logarithmic  backscatter  terms.  The  technique  is,  in 
theory,  not  limited  to  two  dimensional  cases,  but  can  be  used  for 
three  dimensional  grids  as  well. 

Unlike  some  lidar  inversion  processes,  the  tomography  method 
is  not  limited  to  a  fixed  backscatter  to  extinction  relationship. 
This  is  especially  important  for  environmental  conditions  (such  as 
coast  lines)  where  the  aerosols  change  on  short  spatial  and 
temporal  scales . 

All  the  analysis  was  done  for  a  single  flat  scanning  scheme 
where  the  lidar  is  scanned  perpendicular  to  the  average  wind. 
This  scheme  is  achievable  with  many  existing  scanning  lidars,  and 
while  no  effort  was  made  to  find  the  optimal  scanning  scheme  to 
minimize  errors,  wider  scanning  angles  should  reduce  errors. 

The  sensitivity  of  this  algorithm  to  random  movement  restricts 
the  tomographic  applications  of  single  stationary  lidars.  For 
accurate  detailed  analysis,  multiple  spatially  separated  lidars 
will  probably  be  needed.  With  two  scanning  lidars,  the  unknowns 
match  the  number  of  observations  and,  while  the  Inversion  is 
possible,  the  errors  would  be  large,  since  the  errors  for  each 
observation  can  not  be  reduced  by  a  least  squares  fitting  process. 
However,  with  three  lidars,  a  suitable  inversion  should  be 
possible,  since  the  scattering  media  can  be  characterized  from  an 
over-determined  set  of  equations. 
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